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ABSTRACT 

Equilibrium  and  nonequilibrium  molecular  dynamics  (MD)  implemented  on  parallel 
processing  computers  was  used  to  simulate  supercritical  fuel  phenomena  occurring  in  high 
pressure  combustion  devices.  The  coefficients  of  diffusion,  viscosity  and  thermal  conductivity  and 
the  equation  of  state  of  argon,  oxygen,  nitrogen  and  various  alkanes  at  high  pressures  and 
temperatures  were  obtained  via  molecular  dynamics  with  the  obtained  values  agreeing  with  and 
extending  NIST  SUPERTRAPP  code  values.  The  hydrocarbons  ethylene,  butane  and  pentane 
were  simulated.  Vibrational  energy,  dissociation  and  recombination  in  oxygen  and  hydrogen 
diatomic  molecules  were  simulated  with  the  computed  high  temperature  dissociation  rates  in 
agreement  with  published  experimental  measurements. 

TECHNICAL  SUMMARY 

Transport  Properties  and  Equation  of  State.  Molecular  dynamics  has  been  used  to  calculate  the 
coefficients  of  diffusion,  viscosity  and  thermal  conductivity  and  the  equation  of  state  of  pure 
oxygen,  nitrogen,  ethylene,  butane,  pentane  and  hexane  at  high  (supercritical)  pressures  and 
temperatures  .  The  transport  coefficients  are  obtained  by  one  of  several  methods.  Equilibrium 
molecular  dynamics  (EMD)  computes  the  transport  coefficients  by  calculating  the  autocorrelation 
functions  for  a  system  in  equilibrium3.  This  requires  the  computation  of  the  autocorrelation 
functions  over  very  long  time  periods  in  order  to  obtain  accurate  results.  The  shear  viscosity 
coefficient  measures  the  resistance  of  the  fluid  to  a  shearing  force.  For  the  shear  viscosity  the 
integration  is  over  the  non-diagonal  terms  of  the  stress  tensor  and  it  is  a  collective  property  of  the 
fluid.  The  Green-Kubo  formula  is  given  by, 

(i) 

where  Jxy,  an  off-diagonal  term  of  the  stress  tensor,  is  given  by 

*('<)  <2> 

M  ^  1=1  j= 1 

Here  T  and  V  denote  the  temperature  and  volume  of  the  system,  respectively,  and  m  denotes  the 
mass  of  each  particle.  vx,  rij  and  F(rij)  denote  the  x-direction  velocity,  position  vector  and  the 
potential,  respectively,  between  particles  i  and  j.  Statistical  precision  is  improved  by  averaging 
over  all  six  terms  that  result  from  the  stress  tensor: 

V  =  ^  (/Ly  +  ^  +  Vxz  +t*zx+  Myz  +  Vzy)  (3) 

The  thermal  conductivity  coefficient  measures  the  transport  of  heat  in  a  system.  The 
correlation  function  is  obtained  from  the  heat  current  and  it  is  also  a  collective  property  of  the  entire 
fluid.  The  Green-Kubo  formula  is  given  by: 


20000703  020 


0HQ  QUALITY  INSPECTED  4 


2 


1 


(4) 


A  = 


VkbT  0 


](TXM^o+t))dt 


where  T*  denotes  an  arbitrary  component  of  the  heat  current  and  is  given  by 


where  i  denotes  a  unit  tensor.  Again,  statistical  precision  is  improved  by  averaging  over  all  three 
components: 

A  =-(AJ  +Ay  +  AZ)  (6) 

These  functions  are  computed  during  the  simulation.  The  velocity  autocorrelation  is 
computed  outside  the  force  subroutine,  while  the  other  components  of  the  other  autocorrelation 
functions  are  calculated  within  it.  Once  the  correlation  functions  have  been  obtained,  standard 
numerical  integration  techniques  are  used  to  obtain  the  transport  coefficients.  The  statistical 
precision  of  the  self-diffusion  coefficient  is  improved  by  averaging  over  all  the  particles  in  the 
system.  The  viscosity  and  thermal  conductivity  values  are  improved  by  averaging  over  all  the  six 
terms  from  the  stress  tensor,  and  the  three  terms  of  the  heat  current,  respectively.  Figures  1  and  2 
compare  the  MD  calculated  coefficients  of  viscosity  and  thermal  conductivity  for  ethylene  at  10 
MPa  and  a  range  of  temperatures  with  NIST  Database  4  (SUPERTRAPP)  provided  values4. 
Figures  3  and  4  show  the  MD  calculated  coefficients  of  viscosity  and  thermal  conductivity  for 
butane  at  50  MPa  and  Figures  5  and  6  show  the  MD  calculated  coefficients  of  viscosity  and  thermal 
conductivity  for  pentane  at  50  MPa.  The  calculated  values  agree  well  with  the  NIST  values  over 
the  temperature  range  examined  (300-600  K)  for  these  high  supercritical  pressures. 


When  the  calculations  proceeded  to  hexane  the  computation  times  for  the  viscosity  and  the 
thermal  conductivity  autocorrelations  became  excessive.  Nonequilibrium  molecular  dynamics 
(NEMD)  can  be  used  to  calculate  the  transport  coefficients  using  less  computation  time  but  with 
greater  algorithm  complexity5.  NEMD  applies  a  force,  either  diffusive,  momentum  or  energy,  to  a 
system  and  the  rate  of  transport  of  the  desired  quantity  is  computed.  In  order  to  compute  the 
viscosity  we  considered  the  case  of  Couette  flow  in  which  the  fluid  undergoes  sheared  flow  due  to 
boundary  walls  that  are  in  relative  motion.  The  usual  microscopic  definition  of  temperature  in 
terms  of  mean-square  velocity  assumes  that  there  is  no  overall  motion;  any  local  flow  must  be 
subtracted  from  the  velocities  before  using  them  to  evaluate  temperature.  The  same  holds  true  for 
the  velocities  used  in  the  thermostat.  However,  knowing  the  bulk  flow  to  an  accuracy  suitable  for 
use  in  the  equations  of  motion  implies  that  the  problem  has  already  been  solved;  this  circularity  can 
be  removed  by  assuming  the  nature  of  the  flow,  and  only  later  checking  to  see  whether  consistent 
results  are  obtained.  A  less  reliable  alternative  is  to  evaluate  local  flow  by  means  of  coarse-grained 
averaging,  and  then  use  the  results  in  the  equations  of  motion.  Such  an  approach  is  unstable  to  any 
fluctuations  in  the  flow  because  these  variations  are  interpreted  by  the  equations  of  motion  as 
temperature  fluctuations  that  must  be  suppressed. 


In  this  work  we  imposed  the  reasonable  requirement  that  the  MD  flow  obeys  the  linear 
velocity  profile  known  from  the  exact  solution  of  the  continuum  problem.  Assuming  it  is  the  z- 
boundaries  that  are  in  motion,  then  if  the  relative  velocity  of  the  walls  is  yLv  the  shear  rate 
dVj/dzhas  the  constant  value  y.  The  thermostatic  equation  of  motion  is  then, 

r,  =Fi/m  +  a(ri-yrZJx)  (7) 

Where  xis  a  unit  vector  in  the  x-direction,  and  r  =  r(x,y,z)  is  the  position  vector.  The  value  of  the 
Lagrange  multiplier  a  follows  from  the  constant  temperature  constraint: 
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Figure  5.  Shear  viscosity  of  pentane  C5H10,  P  =  50  MPa. 
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a  =  — - - — - = - 

Efc-Fa'X) 
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The  equation  (8)  assumes  that  the  linear  velocity  profile  has  already  been  established.  Creating  the 
initial  sheared  flow  is  most  readily  done  as  part  of  the  initial  conditions,  and  from  the  more  formal 
point  of  view  this  amounts  to  applying  an  impulse  of  the  correct  size  and  direction  to  each  atom  at  t 
=  0.  The  sliding  boundaries,  in  the  form  of  a  special  type  of  boundary  condition,  maintain  the 
constant  shear  rate.  The  constant  temperature  version  of  linear  response  theory  for  this  problem 


provides  an  expression  for  r\  based  on  the  pressure  tensor. 


(9) 


Defining  the  momentum  measured  relative  to  the  local  flow  pjm  =  ri-  yrzi x,  the  first-order 
equations  are  then 

ii=Pi  /m  +  Yrzix  (10) 


Pi=Fi-F2,x+api  (11) 

The  boundaries  are  periodic,  but  of  a  special  form  to  accommodate  the  uniformly  sheared 
flow.  The  idea  is  to  replace  sliding  walls  by  sliding  replica  systems:  layers  of  replicas  that  are 
adjacent  in  the  z-direction  move  with  a  relative  velocity  yLzx,  an  arrangement  designed  to  ensure 

periodicity  at  shear  rate  y.  An  atom  crossing  a  z-boundary  requires  special  treatment  because  the  x- 
components  of  position  and  velocity  are  both  discontinuous  (not  for  the  replica  system  just  entered 
but  relative  to  the  opposite  side  of  the  region  itself  into  which  the  atom  is  actually  inserted).  The 

velocity  change  whenever  a  ±z  boundary  is  crossed  is  +yLzx,  and  the  coordinate  change  is  +dxx, 
where  the  total  relative  displacement  of  the  neighboring  replicas  -  only  meaningful  over  the  range 
(-LJ2.LJ2) -is given  by 

dx=(yLzt  +  LJ2)modLx-L,/2  (12) 

Note  that  because  the  x-coordinate  changes  when  a  z-boundary  is  crossed,  and  additional 
correction  for  periodic  wraparound  in  the  x-direction  may  be  needed.  Interactions  that  occur 
between  atoms  separated  by  the  z-boundary  require  an  offset  value  -dxto  be  included  in  the 
distance  computation. 


To  simulate  the  thermal  conductivity,  a  fictitious  external  field  Fe  of  a  rather  unusual  kind  is 
introduced.  It  has  the  effect  of  driving  atoms  with  a  higher  than  average  energy  in  the  direction  of 
Fe,  while  those  with  a  lower  energy  are  driven  to  the  opposite  direction.  In  other  words  Fe 
generates  heat  flow  and  so,  at  least  for  small  values  of  the  field,  produces  the  effect  of  an  imposed 
temperature  difference. 


The  additional  force  acting  on  each  atom  is  defined  as, 


where, 


e* 


Fe 


The  excess  energy  of  atom  i  over  average, 

The  mean  energy; 

Lennard-Jones  potential  energy; 


ei=\mvl  + 


The  fictitious  external  force; 


(13) 
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f,i  The  force  between  atom  i  and  j  based  on  Lennard- Jones  Potential  model; 

r«  rij  =  ri-rj; 

Na  The  total  number  of  atoms; 

F’  The  additional  force  acting  on  atom  i; 

Here  Ft  has  been  chosen  so  that  in  terms  of  the  heat  current  S  is, 

•F1'  =  VS«Fe 

i 

where  S  is. 


S  =  — 
V 


j  L  i*j 


(14) 

(15) 


N  / 

V  Volume  of  the  system,  V  =  yn, 
n  ’  Number  density,  number  of  atoms  per  unit  volume 
vi  Velocity  of  atom  j 


The  force  conserves  total  momentum  because  ^ Fj  =0,  since  only  relative  distances  occur  in  F-, 

and  assuming  the  force  is  sufficiently  weak  that  the  system  remains  homogeneous,  there  is  nothing 
to  prevent  the  use  of  periodic  boundary  conditions,  exactly  the  motivation  for  devising  methods  of 

this  kind.  If  /  =  Sz,  and  Fe  =  Fe z,  then  the  constant  temperature  version  of  Q  =  lim  lim(/(r))/Fe 
leads  to 

A  =  lim  lim^^  (16) 

F— >o (—>■*>  FeT 

The  thermostat  is  the  usual  one,  based  on  the  total  force,  so  that  the  equations  of  motion  are  simply 

fj  =Fj  +Fi,+«ri  (17) 

Since  the  applied  force  Fe  performs  mechanical  work  on  the  system  the  temperature  rises,  and 
equilibrium  is  never  attained.  To  eliminate  this  problem  a  thermostat  is  included  in  the  dynamics  by 

adding  a  term  apj  to  the  variation  rate  of  momentum  of  every  atom;  constant  kinetic  energy  is 
assured  if  the  value  of  the  Lagrange  multiplier  is 


a  =  - 


F,  +  «F.  4Xfu(ru  •Ft)--i-XfJk(rJk  -F.) 

~j*i  LlSo,  j*k 


•Pi 


(18) 


where  p;  is  the  momentum  of  atom  i,  Pj  =  m\x. 


NEMD  simulations  were  conducted  for  argon  and  ethylene  before  the  grant  period  ended. 
Figure  7  plots  the  thermal  conductivity  of  ethylene  as  a  function  of  temperature  at  5  MPa  pressure. 
Excellent  agreement  between  NEMD  and  NIST  SUPERTRAPP  results  can  be  seen,  especially  at 
the  higher  temperature. 

The  fluid  pressure  for  a  given  temperature  and  density  is  calculated  from  the  virial 
coefficients.  Figure  8  compares  MD  calculated  pressures  for  ethylene  with  experimental  data3. 
Agreement  is  good  except  at  the  highest  density,  which  could  be  improved  by  using  better 
interatomic  potentials. 
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Figure  8.  Ethylene  MD  equation  of  state  calculations  compared  to  experimental  data. 


Vibrational  Energy  and  Dissociation.  As  a  first  step  towards  modeling  chemical  reactions,  the 
vibrational  motions  and  excitation  of  diatomic  oxygen  and  the  subsequent  dissociation  were 
directly  simulated  using  molecular  dynamics6.  The  Morse  potential7  is  used  to  determine  the 
intemuclear  force  between  the  two  atoms  bound  in  the  molecule 


V(r)  =  De~2^r  Tc)-2 De  ^  rc)+E00  (19) 

where  D  is  the  dissociation  energy,  r  is  the  intemuclear  separation,  rc  is  the  mean  bond  length,  and 
Zioo  and  /3  are  constants  computed  from  spectroscopic  data.  Atomic  interactions  outside  the 
molecule  are  modeled  using  the  Lennard-Jones  potential.  For  a  given  temperature  and  density,  the 
MD  calculated  vibrational  energy  distribution  at  equilibrium  demonstrated  excellent  agreement  with 
the  Boltzmann  distribution. 

Dissociation  is  simulated  by  monitoring  the  potential  energy  between  the  two  vibrating 
atoms  in  a  molecule.  When  the  potential  energy  exceeds  a  critical  value  due  to  increased  separation 
between  the  two  atoms,  dissociation  is  judged  to  occur  and  the  Morse  potential  between  the  two 
atoms  is  removed8'10.  Recombination  is  also  possible  by  the  reverse  mechanism.  This  process 
enables  a  system  to  start  in  a  state  of  oxygen  molecules  and  proceed  to  an  equilibrated  state  of 
atomic  and  molecular  oxygen.  Figure  9  shows  the  rate  of  oxygen  dissociation  from  a  single 
unensembled  molecular  dynamics  computation  using  108  molecules  compared  to  that  predicted  by 
kinetic  theory11  using  experimentally  determined  parameters12,13  at  a  temperature  of  6344  K. 
Excellent  agreement  between  the  two  for  both  the  rate  of  dissociation  and  the  final  equilibrium 
value  is  shown. 


Degree  of  Dissociation  T=6344K  v=0.05m3/kg 


Figure  9.  MD  calculated  oxygen  dissociation  rate  compared  to  kinetic  theory  prediction. 
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